#define USE_SSE4 // Nehalem - Sandy Bridge

#include "lamath.h"

#include <iostream>

using namespace lamath;

int
main ()
{

    // Vector Evaluation / Testing
  
    vector v;
    v (0) = 1;
    v (1) = 2;
    v (2) = 3;
    v (3) = 0;

    
    vector w (4, 5, 6, 0);
    
    std::cout << "V1:     \t" << v << std::endl;
    
    std::cout << "V2:     \t" << w << std::endl;
	
    std::cout << "V1 + V2:\t" << (v + w) << std::endl;
    
    std::cout << "V1 - V2:\t" << (v - w) << std::endl;
    
    std::cout << "V1 + 3: \t" << (v + 3) << std::endl;

    std::cout << "V1 - 5: \t" << (v - 5) << std::endl;

    std::cout << "V1 * 10:\t" << (v * 10) << std::endl;

    std::cout << "V1 / 2: \t" << (v / 2) << std::endl;
    
    std::cout << "V1 cross V2:\t" << vector::cross (v, w) << std::endl;
    
    std::cout << "V1 dot V2:\t" << vector::dot (v, w) << std::endl;
    
    std::cout << "V1 Length^2:\t" << v.length_squared () << std::endl;
    
    std::cout << "V1 Length:\t" << v.length () << std::endl;
    
    std::cout << "V1 != V2:\t" << (v != w) << std::endl;
    
    std::cout << "V1 == V1:\t" << (v == v) << std::endl;
    
    std::cout << "+V2      \t" << +w << std::endl;
    
    std::cout << "-V1      \t" << -v << std::endl;
    
    std::cout << "V2 = (5,3,1,9)\t" << (w = vector (5, 3, 1, 9)) << std::endl;
    
    std::cout << "(0,0,1) cross (1, 0 0) = " << vector::cross (vector (0,0,1,0), vector (1,0,0,0)) << std::endl;
    
    std::cout << "\n\n" << std::endl;
    
    // Matrix Evaluation / Testing
    
    matrix m (0.0);
    
    std::cout << "M:" << std::endl << m << std::endl;
    
    for (int i = 0; i < 4; ++i)
    {
	for (int j = 0; j < 4; ++j)
	{
	    m (i, j) = ((4 * i) + j);
	}
    }
    
    std::cout << "M after populating" << std::endl << m << std::endl;
    
    std::cout << std::endl;
    for (int i = 0; i < 4; ++i)
    {
	std::cout << "column " << i << ": " << m.column (i) << std::endl;
    }
    
    std::cout << std::endl;
    for (int i = 0; i < 4; ++i)
    {
	std::cout << "row " << i << ": " << m.row (i) << std::endl;
    }
    
    
    std::cout << "\n-M:\n" << (-m) << std::endl;
    
    matrix t (m.transpose ());
    
    std::cout << "Transpose of M: " << std::endl << t << std::endl;
    
    
    std::cout << "M * Mt:" << std::endl << (m * t) << std::endl;
    
    std::cout << "M == mt: " << (m == t) << std::endl;
    
    std::cout << "M / 4" << std::endl << (m / 4) << std::endl;
    
    std::cout << "Mt * 5" << std::endl << (t * 5) << std::endl;
    
    std::cout << "M * V1: \t" << (m * v) << std::endl;
    
    matrix inv (0.6892133092680031,  2.788855355208904,  2.240785688963923,  4.057479247325732,
		0.43481983566525495, 1.0677941458502471, 3.846393805242511,  1.49122689930281,
		2.9882620097634707,  3.4520774209312766, 2.1713320649601755, 2.945076648179178,
		3.341528216881961,   3.794582616613048,  2.600918779615826,  4.143595773939492);
    
    std::cout << "N: " << std::endl << inv << std::endl;
    
    std::cout << "N^-1:" << std::endl << inv.inverse () << std::endl;
    
    std::cout << "M^-1:\n" << (m.inverse ()) << std::endl;
    return 0;
}
